

 ****************************************************************
 * WELFARE ANALYSIS : FIGURE  2 
 ****************************************************************
 	
	
	use  MxFLS_final , clear		
	  
	 cap drop __*
		  
		   
/* POINT 1 B :  select sub sample   */ 
		
		keep if hhpotentialmigration ==1	 
		keep if statusliving05==1 		
		keep if age >14 
			
/* POINT 2:  indepedent variables   */
		
		#delimit ; 
		global indivlist   female  age age_female   yedu  single   ;  
		global hhlist     		
				ln_socialnonlabincome   						  
				hh_higestgrade3  hh_higestgrade4
				wealthindex  wealthindexsq  
				hhsize     nbage13minus 
				hhrelativeUS hhreturnees
				ARRIVE_age6less    
				LEAVE_age6less 
		; 

		global shocks dead_0205 disease_0205 job_busi_failure_0205 loss_crop_0205 robb_naturaldis_0205  ///
					 dead_0298  disease_0298 job_busi_failure_0298 loss_crop_0298 robb_naturaldis_0298 ;

		global indep_M2     ${indivlist}   ${hhlist}  ${mun_list2} ${shocks}   ;	

		global model_list  M2   ;

		#delimit cr;
	
	*TRIM AT 1% and 99% for all continous variables
	foreach v in cons_attana_0502 { 
      _pctile `v' , p (1 99)
		replace `v' = . if `v' <= r(r1) |  `v' >= r(r2)
	}
	*

 
 * ESTIMATES FOR FIGURE  2
 
    cap drop hhh
	bys folio: g hhh=_n==1

	qui reg   cons_attana_0502  treatment  ${mig_control}    ${hhlist}  nbmaleage3654 nbmaleage1435 nbage55plus   i.id_mun   if rural==1  & hhh==1
	est store m1 
	qui reg hoursagriactivity0502    treatment  ${mig_control}   ${indivlist}   ${hhlist}  i.id_mun    if rural==1  
	est store m0
	suest m1 m0 , cluster(folio) 
	
	forvalues w=0(1)18 {
	qui test [m1_mean]treatment -4*`w'*[m0_mean]treatment =0 ,  matvlc(A)

	global  welfare_gain_`w' =  [m1_mean]treatment -4*`w'*[m0_mean]treatment  
	global  gain_pvalue_`w' = r(p)
	global  gain_std_`w' = sqrt(A[1,1])
	global  diff_cons_gain  4*`w'*[m0_mean]treatment
	
	qui di "marginal productivity  " `w'
	qui di "welfare gain  " ${welfare_gain_`w'}
	qui qui di "standard error  " ${gain_std_`w'}
	qui di "pvalue  " ${gain_pvalue_`w'}
	qui di "diff with pure cons effect   " ${welfare_gain_`w'}
	}
	*
	
	
* FIGURE 2: WELFARE GAIN
	matrix welfare = J(1,4,.)
	forvalues w=0(1)18  {
	
	di "marginal productivity  " `w'
	di "welfare gain  " ${welfare_gain_`w'}
	di "standard error  " ${gain_std_`w'}
	di "pvalue  " ${gain_pvalue_`w'}
	
	di "diff with pure cons effect   " ${welfare_gain_`w'}
	matrix temp= (`w', ${welfare_gain_`w'},${gain_std_`w'}, ${gain_pvalue_`w'}) 
	matrix welfare = welfare\temp
	}
	*	

	cap drop
	svmat welfare
	rename welfare1 mp 
	rename welfare2 esti  
	rename welfare3 std 
	rename welfare4 pval
	g upper_bar= esti+1.96*std 
	g lower_bar= esti-1.96*std 
	
	global xx 
	forvalues w=0(1)15  {
	global xx $xx `w'
	}
	global y -100 -50 0 50 100 150 200 250 300
	
	cap drop a 
	replace mp=2.9 if mp==. 
	g a= 300 if mp==2.9 
	replace a=-100 if mp==3
	replace a=300 if mp==2.9
	
	
	
	twoway (rarea upper_bar lower_bar mp, color(gs13) )  (connected esti mp, color(black) )  ///
			( pci  300 3.1 -100 3.1,  lpattern(shortdash) lcolor(maroon))   ///
			( pci  300 3.15 -100 3.15,  lpattern(dot) lcolor(maroon))   ///
			( pci  300 3.2 -100 3.2,  lpattern(shortdash) lcolor(maroon))   ///
			( pci  300 8.1 -100 8.1,  lpattern(shortdash) lcolor(maroon))   ///
			( pci  300 8.15 -100 8.15,  lpattern(dot) lcolor(maroon))   ///
			( pci  300 8.2 -100 8.2,  lpattern(shortdash) lcolor(maroon))   ///
			( pci  300 13.6 -100 13.6,  lpattern(shortdash) lcolor(maroon))   ///
			( pci  300 13.65 -100 13.65,  lpattern(dot) lcolor(maroon))   ///
			( pci  300 13.7 -100 13.7,  lpattern(shortdash) lcolor(maroon))   /// 
			( pci  0  0  0 15,  lpattern(solide) lcolor(dark))   /// 
			( pci  50 0  50 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  100 0  100 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  150 0  150 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  200 0  200 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  250 0  250 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  300 0 300 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  -100 0 -100 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 
			( pci  -50 0 -50 15,  lpattern(solide) lcolor(dark) lwidth(vvthin)  ) /// 	
			if mp<16, ytitle(Welfare effect) xtitle(Marginal product of self-employed farm labor) ///
			 legend(  si(small) col(2)  symxsize(10)  order( 2 "Point estimate"  3 "Moment of the agricultural" "wage distribution (pesos per hour)" 1 "95% confidence interval"  ) ) ///
			xlab($xx, 	  ) ylab($y, angle(horizontal)  ) 	graphregion(color(white) )	/// 
			 text( -70 3.1 "10th percentile", place(w) color(maroon) ) ///
			 text( -70 13.7 "90th per.", place(e) color(maroon) ) ///
			 text( -70 10.2 "Average", place(w) color(maroon) )
			
	
